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We consider a method for obtaining information on polarization of astronomical objects 
radiation at diffraction limited resolution - differential speckle polarimetry. As an 
observable we propose to use averaged cross spectrum of two short-exposure images 
corresponding to orthogonal polarizations, normalized by averaged power spectrum of 
one of images. Information on polarization can be extracted if object under study can be 
described by model with several parameters. We consider two examples: point-like source 
whose photocenter position depends on orientation of passing polarization and exozodiacal 
dust disc around a star. In first case the difference between photocenter positions can be 
measured with precision of 8 yuas for 2.5-m telescope and 1.2 /ias for 6-m telescope for 
object V = 13 m . For second example method allows detection of discs around central star 
of V = l m with fractional luminosities of 1.8 x 10 -5 and 5.6 x 10~ 6 for 2.5-m and 6-m 
telescope, respectively. 
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Introduction 



Many physical processes give rise to polarization of radiation emitted by a stronomical 



objec ts and measurement of polarization is a powerful observational method ( jTinbergen 



20051 ). In visible light application of polarimetry in combination with high angular 
resolution is promising for study of the following astronomical objects: 1) circumstellar 
environment, 2) Solar System bodies, 3) active galactic nuclei. 

Here we consider a method for obtaining information on polarization of astronomical 
objects radiation at diffraction limited resolution — differential speckle polarimetry (DSP). 
As follows from its title, this method represents a synthesis of speckle interferometry and 
differential polarimetry and presumes the use of instrument combining features of speckle 
interferometer and dual-beam polarimeter. 

Combinat i on of s peckle interferometry and polarimetry has been implemented before 



by iFalcke et al.l (119961 ). Their camera took images in different orientations of polarization 
serially, therefore the data were inevitably affected by differential errors. Nevertheless, 
authors were able to obtain new interesting information on polarization of optical radiation 
of rj Carinae with high angular resolution. 



S chert 1 et al.l (120001 ) constructed dual-beam polarimeter functioning as speckle 
interferometer. A term "speckle polarimetry", which we use here, was coined by authors of 
this work. The instrument was installed on 6-m telescope BTA of Special Astrophysical 
Observatory and worked in NIR band K. Images for each orienta tion of polarization were 
reconstructed separately by bispectrum (ILohmann et al.l . 1 19831 ) and were processed by 



differential polarimetry techniques. These observations should demonstrate that polarized 
flux estimation had a greater precision than total flux estimation, but authors didn't 
indicated it. 

The algorithm of DSP is based on generally accepted model of formation of 
instantaneous image in a focal plane of a telescope: 

G(a) = O(a) <g>T(a), (1) 

where G(a) and O(a) — distributions of light intensity in the focal plane and in the sky, 
respectively, T(a) — instantaneous point spread function (PSF), a — vector of angular 
coordinate in the sky. For convenience let us assume that intensity distribution of object 
is normalized by total flux. 

In Fourier space this equation becomes 

G(f) = 0(f)f(f), (2) 

where G(f) — Fourier spectrum of intensity distribution in the focal plane, O(f) - 
visibility function of the object, T(f) — instantaneous optical transfer function (OTF) of 
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optical system "atmosphere + telescope", / — vector of spatial frequency. OTF fluctuates 
over time due t o perturbation of initially flat wavefront by atmospheric turbulence 
(jGoodmanl Il985l ). As can be seen from the equation, these fluctuations are the source 
of multiplicative noise in G(f), which we will call "atmospheric noise". 



Petrov et al.l (119861 ) have shown that the Fourier spectrum F(f) of image received 



by the detector and normalised by the mean number of photon events K forming it relates 
to the Fourier spectrum of intensity distribution if focal plane in a following way: 



F(f)=G(f)+r)(f), 



(3) 



where rj(f) — spectrum of Poisson noise induced by quantum nature of light, it has circular 
symmetric complex normal distribution with variance equal to K" 1 for all frequencies. 

If the optical system of telescope contains a beam-splitting polarizer, dividing light 
into two orthogonally polarized beams (e.g. Wollaston prism), it will form two images in 
focal plane, corresponding to orthogonal orientations of polarization. In other words, the 
system becomes a dual beam polarimeter. Let us assume that the first beam is horizontally 
polarized (subscript h) and the second is vertically polarized (subscript v). The equation 
([2]) for each of these images becomes: 



G h (f) = O h (f) T h (f), G v (f) = O v (f) T v (f). 



(4) 



After substitution of these equations into (J3]) we obtain: 

F h {f) = (O h (f)f h (f) + Vh (f))e^f\ F v (f) = (O v (f)T v (f) + Vv (f))e^f\ (5) 

Additional phase factors e m ^ h "^' and e m ^ 0v "^ are responsible for displacement of images 
by angles 0y t and 6 V in focal plane, induced by Wollaston prism. Taking account of these 
factors is necessary because in real experiment Fourier spectra of both images will be 
computed in the same coordinate system associated with the detector. 

In speckle interferometry usually a large amount of frames is obtained. After that 
these frames are processed jointly. We also suppose that we have M measurements of 
Fh(f) and F v (f). Let us consider the following value: 



K{f) 



(F h (f)F;(f)) 



M 



(F v (f)F v *(f)) 



M 



K- 1 

v 



(6) 



where (. . . )m means averaging over M measurements. Value (Fh(f)F*(f))M is an averaged 
cross spectrum of imagegj. Expression (F v (f)F*(f)) m in denominator is an averaged 
power spectrum of image and represents biased e stimation of aver aged power spectrum of 
intensity distribution G v (ct). Its bias equals K~ l (jGoodmanl . Il985l ) . therefore denominator 
of (EJ) is unbiased estimation of averaged power spectrum of intensity distribution. For 



1 Petrov et al.l (|1986T ) considered a similar value, but defined for two spectral bands. 
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M > 10 denominator is significantly greater than zero given that the object visibility O v 
is also different from zero. In this case expression fl6]) don't go into infinity. 



Let us assume that the telescope is ideal Th(f) = T v (f) and the object is very bright 
K v <C (F v F*)m- Then, after substitution of (j3J) into (jSj) we get: 



n (f) 



(7) 



One can see that except for the phase factor observable 7Z can be considered as estimation 
of value TZq = Oh/O v , which depends on object properties only. Below we derive more strict 
formula for 7Z and TZq accounting for the instrumental polarization and Poisson noise. We 
also propose a method of estimation of the phase factor. 



Absolute value of TZ was used as an observable by Norris et al ( 120121 ) . They observed 
circumstellar envelopes of red supergiant stars by sparse aperture masking combined with 
dual-beam polarimetry on adaptive optics system NACO/VLT. Success of this work shows 
that parametric analysis of TZ can be quite fruitful. In contrast to that study we consider 
not only the absolute value of TZ , but its argument also, we will call it "phase", by analogy 
with phase of the visibility function. 

We extensively analysed the properties of 7Z, its dependence on instrumental 
polarization of optical system and its noise characteristics assuming that only Poisson 
and atmospheric noises are present. Unfortunately, in general case value TZq doesn't have 
any simple physical meani ng. We propose to use simple parametric model fitting of TZq like 
it was done by Norris et al. (120121 ) . This approach is frequently applied in interferometry for 
analysis of visibility function measure ments. For quantitative evaluations we use images 



simulated using Monte-Carlo method (iMcGlameryl . 1 19761 : lHarding et all Il999l ). 



Theoretical analysis of 1Z properties 



Any real optical system possesses some instrumental polarization. In appendix A we 
consider how it affects the equation of image formation. There we have derived the Fourier 
spectra of intensity distributions in focal plane of dual-beam polarimeter corresponding to 
horizontal and vertical polarization (equations (129]) ) and expression for 1Z (equation (l3Tf ). 

From equation (IHTj) one can see that measured 1Z depends on unknown phase factor 
TZ „ e w((0h-W). This factor can be measured by standard procedure, frequently used in 
differential p olarimetry, — exc hange of images corresponding to horizontal and vertical 
polarization ( iTinbergenl . 120051 ) . This exchange can be performed by half- wave plate, 
installed before Wollaston prism. It allows to rotate the polarization axis of beam by 
it/2. The difference between corresponding measurement 7Z' from original measurement 7Z 
is that phase factor enters it as e m ^ 9v ~ 9h >'f> . Therefore it can be estimated from these two 
measurements by formula 



Jn((O h -0 v )-f) 



vnf)/K'(f). 
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From now forth we assume that the phase factor is known and do not take account of it 
in derivation. Now we consider TZ as estimation of TZq and derive its two basic properties: 
bias and variance. 

Bias of TZ 

In appendix B it is demonstrated for mean of TZ(f) that 

ft(/)=fto(/)(l + Aft(/)). (9) 

This equation says that measurement 71(f) is biased relative to IZo(f). As follows from 
equation f j45l) . the value of bias ATZ(f) depends on polarizing properties of telescope only. 




Fig. 1: Bias ATZ(f) computed using formula ( I45p from appendix B for Cassegrain (upper 
row) and Nasmyth (lower row) foci of the 2.5-m telescope. Left column — real part of 
ATZ(f), right column — imaginary part of ATZ(f). On the axes the normalized spatial 
frequency fX/D is plotted, where A — wavelength, d — diameter of telescope. 

We have computed ATZ(f) for Cassegrain and Nasmyth foci of the 2.5-m telescope, 
results are plotted in Fig. [TJ Real part of ATZ(f) corresponds to amplitude of bias of TZ 
and has typical values of about ±5 x 10~ 4 for both foci. The situation is different for 
imaginary part, which is responsible for the phase of bias of TZ. For Cassegrain focus it 
varies by only ps 10~ 4 , but in case of Nasmyth it shows a significant overall slope with 
amplitude reaching 0.16. This slope is due to inclination of phase in Jones matrix of the 
telescope (see Fig. |H]in appendix A). Therefore we can argue that Cassegrain focus is more 
appropriate for precise measurements of TZq. 

Bias of ATZ(f) can be measured with sufficient precision by means of observation 
of intentionally point-like and unpolarized star given that it and object of interest have 
similar spectra. From now on we assume that we have these measurements and TZ is an 
unbiased estimation of TZq. 
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Variance of 1Z 

It is convenient to decompose the variance of 1Z into variance of its amplitude o~\(f) 
and variance of its phase cr|(/). We estimated values of cr\(f) and cr^(f) using numerical 
simulation described in the next section. We generated 2000 pairs of images corresponding 
to horizontal and vertical polarization. Then for each image we computed Fourier spectrum 
and estimated variances o~\(f) and using formulae from appendix C. 

Let us consider the case of very bright object, in other words let us neglect Poisson 
noise. Results of corresponding computation are given in Fig. |2K as a section of standard 
deviation of cr^if) and cr^(/) along X-direction for Cassegrain and Nasmyth foci. It is 
interesting that variances of amplitude and phase are nearly equal. Also it is noteworthy 
that variances for two foci differ only slightly, what is unexpected taking into account the 
fact that phase of Jones matrix fluctuates over the pupil much more greatly for Nasmyth 
focus than for Cassegrain (see Fig. E] in appendix A). In order to explain this it should 
be recognised that noise of 7Z is largely affected by difference of phases of Jones matrix 
elements Pa and Pjj. Meanwhile for Nasmyth focus most part of phase variation is caused 
by overall slope; this slope leads to bias of 1Z (see previous section), but doesn't affect 
noise. 

In Fig. [2b we give the results of estimations of 0a(/) and cr^(f) accounting for 
the Poisson noise. These values rise at frequencies close to D/X, this is caused by the 
fall of diffraction OTF To(/) at the same frequencies (formula (fTUj) ). The dependence of 
o~A{f) on magnitude can be seen more clearly in Fig. |3] where it is plotted for spatial 
frequency OAD/X. One can see that Poisson noise dominates atmospheric noise for objects 
fainter than V = — l m , i.e. for almost all astr onomical objects. I n domain of Poisson 



noise prevalence the approximate formula from (IPetrov et al.l . 119861 ) describes behaviour 



of variance reasonably well (thin dashed line in Fig. [3]): 

oi(f)=(M\T(f)\ i f {f)KL- i y\ (10) 

where /(/) — visibility function of object (here we assumed that object appearance doesn't 
depend on polarization significantly), L = 2.3(D/ro) 2 — mean number of speckles in the 
image, r — Fried parameter. In domain of Poisson noise prevalence difference between 
noises for Cassegrain and Nasmyth foci essentially disappears, consequently hereafter we 
consider the Nasmyth focus only. 

Monte-Carlo simulation of DSP 



In order to estimate the expected perf ormance of the DSP we simulated the pro cess of 
measurements using Monte-Carlo method (iMcGlameryl . Il976t lHarding et al.l . Il 9991 ) . This 
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Fig. 2: Section of a a (solid lines) and (dashed lines), on the horizontal axis the 
normalized spatial frequency fX/D is plotted, a — computation without atmospheric and 
with Poisson noises, black and grey lines — for Cassegrain and Nasmyth foci, respectively; 
b — computation with atmospheric and Poisson noises and for different magnitudes of 
object, values are given in figure. Computation details are given in text. 

method is based on generation of a large number of instantaneous images, disturbed by 
the atmosphere, followed by simulation of subsequent processing. 

We modelled the atmosphere as a set of infinitely thin turbulent layers. Each of these 
layers is defined by altitude, turbulence intensity, wind speed and direction. For each layer 
we generated a realisation of random phase screerj§ constituting N x N matrix with von 
Karman power spectrum: 

F,(f) = 0.0229 r~ 5/3 (/ 2 + Lfy ll/ \ (11) 

/ — absolute value of spatial frequency vector /, r — Fried parameter, corresponding to 
turbulence intensity of this layer, Lq — outer scale of turbulence, we adopted typical value 
of Lq = 25 m. 

2 Phase screen is an idealized infinitely thin optical component, which affects only phase of passing 
wavefront. 
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Fig. 3: Dependence of a a on magnitude V for / = OAD/X, Nasmyth (thick solid line) 
and Cassegrain (thick dashed line) foci. The leftmost points are computed for very bright 
object. Computation details are given in the text. Thin dashed line stands for approximate 
estimation of Poisson noise computed using formula (ITU]) . 



Originally flat wavefront propagated these phase screens, starting with the highest, 
in output we obtained disturbed wavefront. Then we took into account that instrument 
is dual-beam polarimeter. In order to d o this we made t wo co pies of wavefront. Then in 



accordance with Jones matrix definition (IBorn and WolJ . 1 1 9 T3h we multiplied one copy by 
Pa and another by Pp, see appendix A. 

Calculation of Pa and Pp was made by use of Zemax software for three optical 
systems: Cassegrain and Nasmyth foci of the future 2.5-m telescope of SAI and for 
primary focus of the BTA. For the 2.5-m telescope the central obscuration is 0.43, for BTA 
this parameter equals 0.33. We assumed that mirrors are coated with bare aluminiunfl 
Amplitude and phase of calculated Jones matrix for the 2.5-m telescope are presented in 
Fig. 0in appendix A. 

After this we computed disturbed images and added the Poisson noise to them. 
We have taken into co n sidera tion two-fold multiplicative noise of EMCCD detector 



( iHynecek and Nishiwakil . 120031 ). optimal for speckle interferometry. Thus we have 
accounted for all effects of interest for our problem: atmospheric disturbance, instrumental 
polarization and Poisson noise. 



3 Transparent protective coating of mirrors can affect the instrumental polarization greatly. Therefore 
for mirrors of real telescope instrumental polarization can be substantially different from our calculations. 
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Let us discuss some parameters of the simulation. Computation was performed for 
V band, some parameters were also computed for / band. Total efficiency in V and / was 
0.54 and 0.43, respectively (with consideration of transmission of atmosphere, telescope 
and instrument optics, quantum efficiency of detector). We adopted value of 30 ms for 
exposure, what is typical for speckle interferometric observations. For objects fainter than 
l m we used minimal period of exposures 30 ms, as long as in this domain uncorrelated 
Poisson noise dominates. For bright objects period was increased to 120 ms for adequate 
account for atmospheric noise. 

Uncorrelatedness of Poisson noise allowed us to compute quantitative performance 
characteristics — variances of estimations — for small number of frames, e.g. 100, and 
then adapt them for long series by multiplication by square root of frame number ratio. 
Standard series duration was adopted as 1 hour, what is much longer than typical series 
in speckle interferometry. For correct processing of speckle interferometric measurements 
atmospheric conditions (Speckle Transfer Function — STF) should be constant throughout 
the series. For DSP there is no such requirement. 

As a model of atmosphere we adopted two turbulent layers of equal intensity, yielding 
seeing 0.91" in V band . This is median seein g for Mt. Shatdzhatmaz, supposed location 



of the 2.5-m telescope (IKornilov et al.l . 120101 ) . Both layers move with speed of 10 m/s in 



opposite directions, what provides maximum temporal uncorrelatedness of atmospheric 
noise. This rough model is sufficient as long as we are not interested in temporal or 
isoplanatic characteristics of image. For convenience of comparison computation for BTA 
was performed for the same model. 

Resulting series of images and F v were used for estimation of 1Z by formula (jSJ). 
Its variance was estimated by formulae from appendix C. 

Parametric analysis of 1Z 

Fitting measurements of 1Z with a model of a few parameters is the simplest way of 
extraction of information from it. We consider a special case of such a model, having one 
parameter: 

$(/,p) = $o(/)+P$i(/). (12) 

In this case the problem resolves itself into finding of parameter p, for which function 
$(/,£>) describes the observations in the best way. We estimate the parameter by 
minimization of sum of squared deviations of model from measurements, taken with weights 
inversely proportional to estimated variances of measurements. 

Error of p estimation, or equivalently its expected variance a^, defin es the DSP 



efficiency. For given linear model it is derived by the following expression ( iKuzmenkov 
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(119851 ) considered more general case): 



a; — at 



Summation is being performed over {f c /fo) 2 independent measurements, where f c = D/X 
- cutoff frequency, of — variances of measurements, they can be estimated using formulae 
( j47l) and ( 1481) . f i — coordinates of measurements of 1Z in Fourier space. 



Now we consider two parametric models: 1) point-like object, whose photocenter 
position depends on orientation of polarization; 2) point-like object and faint polarized 
envelope. 

Point-like object, whose photocenter position depends on orientation of polarization 

Suppose that we have an object whose appearance depends on orientation of 
polarization Oh ^ O v , therefore TZ ^ 1. Let us denote typical angular extent of this object 
as 7, what corresponds to spatial frequency domain / 7 « I/7. For frequencies |/| <C / 7 
amplitude and phase of TZq can be represented by Tailor series with small parameter /// 7 . 
In this domain variation of amplitude is dominated by quadratic term and variation of 
phase — by linear term. If object demonstrates high asymmetry of polarized and/or total 
flux then TZq phase changes significantly at frequencies ps f 7 . In this case in the domain 
I /| <C / 7 variation of argument is much larger than variation of amplitude, moreover, most 
part of this variation consists in slope. Consequently, for 1Z in domain |/| <C / 7 one can 
write: 

n ^exp{nr(0f)}, (14) 

where 6 — vector of difference between photocenters of images corresponding to different 
orientations of polarization. In reality frequencies |/| <C / 7 are sampled when the object is 
much smaller than diffraction resolution of telescope. Thus becomes the most robustly 
defined parameter of object. 

For the described type of object it is convenient to model the phase of TZ. In this 
case the model is <&(f,6 x ,0 y ) = iv(f x 6 x + f y y ), where 9 x ,6 y — components of vector of 
photocenter position difference. We take the 6 X as a sought quantity, 9 y is assumed to be 
known. This simplification is valid as long as mutual estimation of 9 X and 9 y is uncorrelated. 
In terms of the linear model (fl2"j) 3>i(/) = itfx- Taking this into consideration we can 
rewrite equation (|T3l) for the parameter 9 X . 

^£^#. (15) 

where f xi — X-component of vector f { , phase noise a'L = cr^/J is derived using method 
from appendix C. 
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Using this expression we computed expected variance of angle 6 X for adopted 
conditions of simulation and plotted it in Fig. 0K It can be seen that expected variance 
doesn't de pend on seeing. Th is feature of differential speckle interferometry was noted 



earlier by (IPetrov et all 119861 ) . 



The figure shows that for object V = 13 m , precision of estimation of difference 
between photocenters corresponding to horizontal and vertical polarizations is 8 //as for 
the 2.5-m telescope and 1.2 /xas for the 6-m telescope. Such precision significantly exceeds 
typical precision of differential narrow-field astrometry at 5-8 m telescope s, which reaches 



100 / x as and 150-200 /xas with and without adapti ve optics, respectively (ICameron et al. 
20091 ; IPravdo and Shaklanl . Il996l : lLazorenkol . l2006h . 



The detection and measurement of differential astrometric signal can be particularly 
productive for the study of BL Lac objects. These active galactic nuclei (AGN) periodically 
show increase in polarization degree of 10-20% on timescales of 50-100 days. It would 
be interesting to measure the polarized flux photocenter movement accompanying these 
events. 



Exozodiacal disc 



Various processes associated with stellar evolution frequently give rise to the 
occurrence of circumstellar matter, e.g. dust disc. Radiation of dust has two main 
components. The first is thermal IR radiation of dust heated by nearby star. At the moment 
observations and analysis of IR radiation of dust is the main source of information about its 
distribution in circumstellar space, temperature and chemical composition. This is caused 
by the fact th a t in I R the contrast between disc and star is not so large as in visible 
(jdi Folco et all 120071 1 . The second component of dust radiation is scattered visible light 
of star. Imaging of dust discs in scattered light, especially in multiple photometric bands 



is a p owerful tool of dust diagnosis, which supports IR observations ( lAbsil and Mawet 



20101 ) . Main difficulty of circumstellar dust observations is high contrast between star 
and envelope — 10 -5 and higher. However scattered light is polarized what simplifies the 
problem of its detection and characterization. In this subsection we determine limiting 
contrast between star and disc for DSP using numerical simulation. 

Let us assume the following simple model of object: exozodiacal disc, analogous to 
solar one, but with total luminosity z times higher, rotating around r Ceti (L = O.52L , 
distance 3.65 pc, magnitude V = 3.5 m ), visible from the pole. Parameter z is a measure 
of disc's "m agnitude" relat i ve to solar one, a term "zodi" was adopted in l iterature as a 



unit for it (IRoberge et all 120121 ) . The disc was simulated with ZODIPIC (IMoran et al. 



20041) package for I DL. This program makes use of data on the Solar zodia cal disc from 
( Moran et al. . 2004 ). Scattering by dust is modeled in accordance with paper (Hong], 1985 ). 
We augmented the code with the computation of Stokes vector components Qd(oc) and 
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Fig. 4: a — dependence of expected differential astrometry precision on magnitude V. b 
- dependence of expected precision of z estimation on magnitude V for exozodiacal disc 
analogous to disc around r Ceti. z — flux of disc relative to flux of solar zodiacal disc. 
Along right vertical axis the contrast k of disc is plotted. For both figures: solid lines - 
r =12.5 cm (median conditions at Mt. Shatdzhatmaz), dashed lines — r = 20.5 (10% of 
best conditions at Mt. Shatdzhatmaz). Black lines — 2.5-m telescope, grey lines — 6-m 
telescope. Dotted line stands for median conditions, 2.5-m telescope and I band. 



Ud(oi), corresponding equations also were taken from ( iHongj . Il985h . Obtained distributions 
of Stokes parameters I, Q, U are given in Fig. |5j 



Visibilities Oh and O v for this type of objects are 

Oh = Istar + z Id + z Qd, O v = I star + Zlj, — zQ d , 



(16) 



where I star — visibility function of the star, Id and Qd — Fourier spectra of Stokes vector 
distribution of the disc. The Qd for our model of disc is also plotted in Fig. In spite of the 
fact that disc angular size is formally below diffraction limit of the 2.5-m telescope, evidence 
of its presence in Qd extends into Fourier space domain available for measurements with 
the 2.5-m telescope (i.e. at frequencies smaller than cutoff frequency). The possibility to 
obtain information about s ource features smal l er tha n diffraction resolution of telescope is 
known as super-resolution ( jMarti-Vidal et al.l . l2012h . 
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Value IZq for this model becomes 

^, _ I star + Zl d + zQ d 

I-Cq - ~ ~ [I ( ) 

istar + Zld — ZQd 

Let us assume that total luminosity of envelope is smaller than stellar luminosity 

Id ^ Istar'- _ _ _ 

Wo «(i + i^)(i + ^)«i + ?^. as) 

V Istar J V Istar J I star 

Now we take into account that star is point-like source and its visibility I s t ar — P 

TZ ^1 + 2zQ d . (19) 

Thus parametric model, approximating observed TZ, is: 

<f>(f,z) = l + 2zQ d (f). (20) 

The aim of fitting is finding of z. Another consequence from this equation is that it is 
possible to obtain Fourier spectrum of polarized flux of circumstellar envelope and then 
its image. However image has a very large number of unknowns (pixels intensities), each 
of them will be estimated with low precision. That is why for faint discs model with small 
number of parameters (in our case just one) is more appropriate. 

In this case we will be interested in amplitude of 1Z only, because the source is 
centrosymmetrical, therefore phase of 71 equals zero. Equation ( [13"]) for the error of z will 
be: 

1 V- ±Ql 



al ^ a 2 - ' 



(21) 

'At 

where Qdi = Qd(fi), noise of amplitude a 2 Ai = o\(fj) is being estimated by formula (l47j) . 

Parameter z estimated from measurements is a normally distributed random number, 
therefore probability of detection of disc z = a z is w 70%. Although this threshold 
corresponds to certain model, it provides a glimpse of method performance in observations 
of circumstellar environment. 

Estimations of detection limit for the 2.5-m and the 6-m telescopes and different 
seeing conditions are presented in Fig. Hb. One can see that in this case seeing also weakly 
affects performance. 

At the 2.5-m telescope DSP will allow detection of discs z = 700 zodi around V = l m 
stars, at 6-m telescope — z = 220 zodi (absolute contrast between the disc and the star 
1.8 x 10~ 5 and 5.7 x 10~ 6 , respectively). Exozodiacal disc around Vega is one of brightest 
discs, in IR its luminosity reaches 3000 zodi. In visible luminosity of the disc is likely to 





Fig. 5: The model of exozodiacal disc around r Ceti, computed using ZODIPIC (for details 
see the text). Upper row: Stokes parameters U and Q, lower row, left panel — total intensity 
/, circles diameters correspond to 1.22X/D for D = 2.5 m (large circle) and D = Qm (small 
circle). Lower right panel — Fourier spectrum Q^, circles stand for domains available for 
measurements with the 2.5-m (small circle) and the 6-m (large circle) telescope. 



be lower ( Absil et al. . 2006 1^1. Thus we can conclude that brightest exozodiacal discs are 
available for study with DSP at the 2.5-m and the 6-m telescopes. 

Exozodiacal discs are only special case of circumstellar envelopes, they are typical for 
relatively old main sequence stars. In other situations contrast between star and envelope 
is much lower. E.g. for envelope observed by (INorris et all 120121 ). fraction of scattered light 
~ 10%, therefore it can be de tected by DSP ea sily. Young stars also frequently possess 
strong dust envelopes (see e.g. (IKrist et all 120051 ) ). These objects are of great interest and 
results of this subsection can be used for estimation of possibility of their observation with 
DSP. 



Practical aspects of DSP 

Specifics of DSP impose some constraints on device implementing this method. Let 
us discuss the most important of them. 

4 The data on exozodiacal discs luminosity in visible are quite scarce due to mentioned difficulties of 
high contrast. 
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Fourier spectra used for calculation of TZ can be correctly estimated if pixel angular 
size is less than X/2D. EMCCD detector is preferable as long as it has negligibly small 
readout noise and allows to make short exposures at sufficiently high frequency. 

Theoretical analysis shows that TZ is biased estimation of value TZo = Oh/O v , which 
depends on object properties only. Bias of TZ depends on polarizing properties of telescope 
and instrument, for measurement of it we propose two stages of calibration. 

The first stage is exchange of images corresponding to horizontal and vertical 
polarizations, what can be implemented with half-wave plate rotating polarization axis 
by 7r/2. This procedure will allow to precisely measure phase factor from equation (J7J), or 
equally the difference between deviation angles of Wollaston prism. Differential aberrations 
occurring in optical system after half-wave plate also will be contained in this phase factor. 
Thus it makes sense to install half-wave plate as early as possible in optical system. 

The second stage of calibration is required for removal of differential aberrations 
induced by optical elements before half- wave plate (e.g. telescope mirrors). This calibration 
consists in measurement of TZ of some reference star. This star should have similar spectrum 
and be brighter than object of interest by 4 — 5 m (when it is possible), otherwise calibration 
measurement will significantly increase the error of final estimation. Appearance of 
reference star shouldn't depend on orientation of polarization axis. This condition will 
be fulfilled more than enough for single main sequence star of spectral class later than AO 
having angular size less than « 10 mas. 

Differential aberrations give rise to not only bias but atmospheric noise as well. 
Therefore they should be minimized in design phase. This can be achieved by ensuring 
that paths of beams corresponding to horizontal and vertical polarizations are as close 
as possible to each other. In order to do this Wollaston prism should be installed as last 
optical element in system. In this case the most part of differential aberrations will be 
induced by the prism itself due to birefringency of its material. 

Deviation of any surface of the prism intersected by beams from the plane by A 
results in differential phase aberration k(n Q — n e )A, where n Q and n e — refractive indices 
corresponding to ordinary and extraordinary beams. Value n — n e for calcite is —0.172, 
for quartz +0.009. Thus the use of quartz prism is more appropriate, as long as surface 
quality requirements are relaxed for quartz in comparison with calcite. However deviation 
angle for quartz prism is significantly smaller. 

Let us estimate acceptable level of differential aberrations. Differential astigmatism 
(similar to plotted in Fig. O in appendix A, lower right, phases of Pa and Pp) Z 5 = 
0.011, what corresponds to amplitude of phase variation 0.054 rad, results in increase of 
atmospheric noise by factor of 3 in Cassegrain focus. Corresponding deviation of prism 
surface from plane is A/20 for calcite and A for quartz. In this case atmospheric noise 
mainly affects observations of bright objects V < 4 m , because for them Poisson noise is 
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sufficiently small (see Fig. [3]). If we restrict ourselves to faint objects, requirements for 
differential aberrations are substantially relaxed. 

In analysis being performed in this work we take into account three effects: 
instrumental polarization of telescope, atmospheric disturbances and finiteness of number 
of detected photons. Measurements made with real instrument will probably be affected by 
other systematic errors induced by other factors. However some of these errors will likely 
to be eliminated by described calibrations. 



Authors of (jCanovas et al.l . 120111 ) achieved significant success in removing of 
systematic errors in differential polarimetry with Extreme Polarimeter (ExPo) instrument 
at 4.2- m WHT telescope. This instrument is quite similar to one supposed by us, it obtains 
a large series of short-exposure images taken in orthogonal polarizations simultaneously. 
ExPo exchanges the images after every exposure for elimination of systematical errors. It 
is probable that this approach will be effective for DSP also. 

Conclusions 



Differential polarimetric measurements with high angular resolution are much more 
precise than absolute what provides additional opportunities of study of astronomical 
objects. Improved precision arises out of the fact that many factors disturbing wavefront 
affect its polarization components equally. Atmospheric turbulence limiting angular 
resolution of ground based telescopes is one of these factors. DSP allows to significantly 
reduce atmospheric noise influence and obtain information about polarization properties 
of object's radiation with diffraction resolution. 

As an observable for DSP we have adopted value 1Z — averaged cross spectrum 
of two short-exposure images corresponding to perpendicular orientations of polarization, 
normalized by averaged power spectrum of one of images. In order to analyse the properties 
of 1Z noise we perform numerical simulation of DSP using Monte-Carlo method. The main 
outcome of this simulation is that Poisson noise dominates atmospheric noise for stars 
fainter than V = — l m given that there isn't any differential aberrations apart from induced 
by the telescope. 

Measurements of 1Z can be interpreted by parametric model fitting. Using numerical 
simulation we study the performance of method for two examples of such model. In first 
example we consider point-like object whose photocenter position depends on orientation 
of polarization. In this case difference between photocenter positions for object with V = 
13 m can be measured with precision of 8 /ias and 1.2 /ias at the 2.5-m and the 6-m 
telescopes, respectively. Such astrometry precision would provide interesting opportunities 
of investigation of some astronomical objects, e.g. AGN. In second example we evaluate 
capabilities of detection and characterization of circumstellar dust discs using idealized 
model of such disc. It is found out that for V = l m star discs with luminosity of 1.8 x 10~ 5 
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and 5.6 x 10~ 6 relative to the star can be detected at the 2.5-m and the 6-m telescope, 
respectively. Such contrast is typical for brightest exozodiacal discs and dust envelopes 
around young stars. Performance of DSP depends on the seeing weakly. 

It is interesting to compare capabilities of DSP with other existing and planned 
polarimeters. At the moment one of the pressing observational problems in astronomy is 
the detection of faint components (mainly exoplanets) in close vicinity of bright stars. 
There are several existing; and planned i nstruments de signed to so l ve th is problem: 
HiCIAO /SUB ARU jHashimoto et al.L boilh. NACO/VLT dNorris et alll2012h. S PHERE 



IRDIS/VLT (Dohlenetal 



2008h. SPHERE-ZIMPOL/V LT flThalmann et al 



2008|) 



GPI/Gemini (Macintosh et all . 120061 ). MMT-POL/MMT flPackham et all . 120101 ). Some 
of them will be able to reach contrast of 10~ 8 . All these instruments are very complex and 
expensive devices having differential polarimetry mode and fed through adaptive optics. 
Meanwhile all of them except ZIMPOL are designed for NIR. 

ExPo instrument works in visible, but its processing algorithm operates with images, 
therefore its angular resolution 0.5" is limited by the seeing. Angular size of pixel of this 
camera is larger than X/2D, therefore the data obtained with it cannot be processed with 
DSP. 

Thus it can be seen that DSP has its niche — relatively high contrast even in visible 
at small angular separations from the star and high precision of differential astrometry. 

Author is grateful to Tokovinin A. A. for discussion of speckle interferometry issues. 
Comments from Kornilov V.G. and anonymous referee helped to substantially improve the 
paper. Author acknowledges support from "Dynasty" foundation. 

Appendix A. Polarization in focal plane of telescope 



Vector Sf(h,v)(f) of Fourier spectra of Stokes parameters distribution in the focal 
plane of the telescope re l ates to similar vector S(f) defined in the sky in the following 
way (Almeida and Pilletl . 1 19921 ) : 



V(M) 



(f) = M f(h , v) (f)S(f). 



(22) 



This vector equation is written down for both images corresponding to horizontal (h) 
and vertical (v) polarizations. Mfr^if) are Fourier spectra of Mueller matrices of the 
telescope. Equation (I2"2"j) is valid as long as Mf^,v){f) doesn't depend on direction in the 
sky. This condition may be considered fulfilled because we deal with narrow fields ~ 2". 
Equation (122]) constitutes some generalization of the common equation of image formation 
in Fourier space, in this sense Mf^,v){f) is generalized OTF. 



The method of computation of matrix Mf^,v){f) is given in (lAlmeida and Pillet 



19921 ) and (lAzzam and Basharal . 119871 ). Here we briefly reproduce it: 



M /(M (/) = BT h>v (f)B- 



(23) 
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where matrix B is defined (lAzzam and Basharal . 119871 ) as 



B 



(l 
1 
1 

\0 i 



1\ 

-1 

1 
-% J 



(24) 



Matrices Th(f) and T v (f) also have size 4x4. Given that Wollaston prism 
is ideal all elements of T/j(/) equal zero except for the elements of the first row: 
TAA*(f),TAB*(f),TBA*(f),TBB*(f)- The same is valid for matrix T v (f), however in this 
case the last row elements are non-zero: Tcc*{f),TcD*(f),TDC*(f),TDD*(f)- Elements of 
these matrices can be computed as follows: 



T XY *(f) = Tl- 1 / P x (x)P*(x + Xf)exp{-i(<f>(x)-<f>{x + \f))}dx, 



(25) 



where X, Y run over A, B, C, D, <f>(x) — variations of phase induced by th e atmosphere, II 
total area of aperture. For polarization computations a Jones matrix (IBorn and Wolf . 



19731 ) constitutes analogue of the aperture function: 



P(x) 



Pa(x) 
Pc{x) 



P B (x 
Pd(x 



(26) 




We use elements of this matrix in formula (1251) . The Jones matrix defines the relation 
between Jones vectors before (J) and after (J 7 ) optical element: 



(27) 



Outside the aperture all elements of this matrix equal zero. 

Amplitudes and phases of the Jones matrix for Nasmyth focus of the 2.5-m telescope 
are plotted in Fig. EJ The computation has been performed with Zemax software. 
Evaluation shows that the amplitude of diagonal elements Pa, Pd has mean over the 
pupil Rj 0.87 and fluctuation m 0.2 — 0.3%. In contrast to them non-diagonal elements 
have small mean 1% and fluctuations of the same order 1%. Thus the amplitude is much 
smaller for non-diagonal elements then for diagonal. Phases of the diagonal elements differs 
by w 0.1 rad. Results for the Cassegrain focus are also given in Fig. El Absence of oblique 
reflection gives rise to much smaller variations of the P matrix elements. For diagonal 
elements Pa, Pd variations of the amplitude amount to 0.05%, for non-diagonal — 0.2%. 
Difference in phases of diagonal elements also is much smaller — 0.01 rad. 



Inserting all elements of f l2"5"j) into equation ( 12"3"|) and using expression (I2"2"j) we obtain 
Fourier spectrum of Stokes parameters distribution in the focal plane Sf(h,v)(f)- Taking 
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amplitude 



phase 



■0.8775 



0.013 



t/a 
o 



2 




0.8575 




0.13 



0.08 



«3 
O 

.2 



bo 

t/3 

u 



W 



'0.8935 
10.0042 



0.0042 



J 




3.14 







i m i 

3.0045 — '-3.14 

S.14 ^■^—^^■0.0045 

I 

0.0045 



Fig. 6: Jones matrices. First and second row — computation for the Nasmyth focus, third 
and fourth — for the Cassegrain focus. First and second column — amplitudes (from left to 
right, from top to bottom: Pa, Pb, Pc, Pd), third and fourth — phases (from left to right, 
from top to bottom: Pa, Pb, Pc, Pd)- 

into account that we measure the intensity only, we obtain resulting expressions for its 
Fourier spectrum: 



G h = (T AA * + T BB *)I + (Taa* - T BB *)Q + (Tab* + T BA .)U + i(T BA * - T AB *)V, 
G v = (f cc * + Tdd*)T+ (fee* - f DD *)Q + (fen* + f DC *)U + l(f DC * ~ f C D*)V. 



(2* 



Here the /, Q, U , V — Fourier spectra of Stokes parameters distribution of the object. As 
long as evaluation shows that non-diagonal elements Pb, Pc of the Jones matrix have first 
order of vanishing relatively to diagonal Pa, Pd, equations (128|) can be simplified given 
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that fraction of polarized flux is small / ^> Q,U,V: 

G fc (/)«r ilA .(/)(7(/) + g(/)), 

G v (f)*T DD 4f)(I(f)-Q(f)). 

After substitution of these equations into ([3]) and taking into account unknown phase 
factors we obtain Fourier spectra of detected images 



F h (f) = [T AA *(f)(I(f) + Q(f)) + r]h (f)]e i ^f\ 
F v (f) = [TW (/)(!(/) - Q(f)) + r j v (f)]e i ^\ 

Substitution of (130]) into ([6]) gives in turn: 



(30) 



[Taa* (I + Q) + Vh] [Tdd* (I-Q) + w] * ) e M(« h -«-)-/) 

K=^-z 3—3 3—3 ^ , (31) 

[T DD . (I — Q) + Vv ] \Tdd* (I - Q) + Vv] ) ~ K-i 



M 

where dependence on / is omitted for brevity. 

Appendix B. Evaluation of 1Z mean 

The value 1Z represents ratio of sample means computed from M measurements. From 
statistics it is known that sample mean of some value is a normally distributed random 
number with mean equal to mean of this value. Therefore, given that M is sufficiently 
large, TZ is also normally distributed random number with mean equal to the ratio of 
means of values FhF* and F V F* — K~ x . Let us denote this mean as TZ and let us estimate 
it. The equation (131]) is taken as a starting point. 

Let us introduce the following notations: 

NM --7(lffkff (32) 

where /(/) + Q(f) and 1(f) — Q(f) are constant in time. 

Given that fluctuations of Pa{x) and Pd{x) are small is comparison with the values 
themselves the Taa* and Tdd* can be transformed (we perform derivations for Pa{x), for 
Pd(x) they are analogous): 

P A (x) = P (x) + AP A (x), (33) 

where Pq(x) is the aperture function in the usual sense, it equals one inside pupil and zero 
outside, APa{x) ^ 1. Inserting this into integrand of (125]) we get (neglecting the cross 
term, because it has second order of vanishing): 

Taa-W = r atm (/) + ATaa* (/)) (34) 
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where T a ^ m (/) is OTF of system telescope + atmosphere in the usual sense: 

?atm(/) = n- 1 J Po(x)P*(x + Xf) exp {-i(<f>(x) - <f>{x + A/))}dac. (35) 
AT A A*{f) is some correction factor similar to OTF: 

AT AA *{f) = IT 1 J P (x)P*(x + Xf) [AP A (x) + AP* A {x + Xf)] 

x exp {-i{<j){x) -<j>(x + Xf))}dx. ( 36 ) 

II is a total aperture area. Substituting (134")) into (131]) and after development we obtain 



T atm + AT AA * + JV fc [T atm + ATdd* + N v ] 

K = -r-^ , (37) 

[T atm + ATdd* + N v ] [T atm + ATdd* + A^] \ - tf-i 



where dependence on / is omitted for brevity. Hereafter angle parenthesis means averaging 
over assembly. 

Now we are going to recast (l37p : remove parenthesis, take into account that 
ATaa*, ATdd* <C T a ^ m , (Nh) — (N v ) = (as long as they are circularly symmetric 
complex normal random numbers) and (T a ^ m ) = (ATaa*) = (ATdd*) = for frequencies 
f > vq/X (absence of high spatial frequencies in long exposure image): 



K(/) = ftofi+ ( 5 tm . Ar) 

( T atm T a t m / - 



(3* 



where AT = ATaa* — ATdd*- Let us develop it explicitly: 

AT (/) = IT 1 J P (x)P*(x + Xf)AP(x, f) exp {-i{<f>{x) - 4>{x + Xf))}dx, (39) 
where AP: 

AP(x, f) = APa(x) + AP*> + A/) - AP D (aj) - AP* (a; + Xf). (40) 
For numerator of the second term in fl38l) one can write: 



(T* tm AT) =n- 2 / / P (x)P*(x + Xf)AP(x, f) 

J J (41) 

x (exp {-i{(f>{x) - <f)(x + Xf) - <f>(x') + <j){x' + Xf))})dxdx' . 



Averaged expression in this integral equals a5(x — x') for frequencies / ^> r /A, 
where a is at mospheric coherence area, which equals 0.342rQ for the Kolmogorov turbulence 



(iKorfi 119731 ). Using this fact we obtain: 

(f * tm AT) = (ZFI)- 1 J P (x)P*(x + Xf)AP(x, f)dx, (42) 
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where L — total number of speckles in the image. 



On the same assumptions one can show that (jKorfj, Il973l ) 



( T atm^atm> = (^n)" 1 J Po(x)P^x + A/)das, 



(43) 



where integral is OTF of the system in the absence of the atmosphere T (f). Ultimately 
forft(/): _ 

K(f)=Ho(f)(l + AK(f)), (44) 

where 

AK(f) = J P (x)P*(x + Xf)AP(x, f)dx I J P (x)P*(x + Xf)dx. (45) 

Thus ATZ(f) is AP(x, /) averaged over the area defined by Pq(x)Pq(x + A/). 

Appendix C. Evaluation of 1Z variance 

Here the method of evaluation of variance of TZ from simulated data is given. As can 
be seen from (J6]), TZ can be considered as a function of three random numbers X, Y, Z: 

where X = Re(F h F*) M , Y = Im(F h F;) M , Z = {F V F*) M . 

For further use it is convenient to decompose the variance of the TZ into variance o\ 
of its amplitude and variance cr| of its phase. On the assumption that standard deviations 
of X, Y, Z are small relatively to the values themselves, one can express a\ and cr^ as 
follows: 

2 2 X 2 2 Y 2 2 w2 2 2XY 2 2X 2 2Y 

°A - °XX Z ,2 W 2 + a YY Z ,2 W 2 + °>Z^4 + G XY Z ,2 W 2 ~~ a XZ^ ~ °YZTfii ( 4 <) 

2 &xx y2 + Q'yyX 2 - 2o~XY XY /a a\ 

a t> = m ' (48) 

where crxx-> a XYi e ^ c - are elements of the covariance matrix of values X, Y, Z. Here we 
use delta method and introduce two auxiliary values for convenience: Z' = Z — K~ x and 
W = VX 2 + Y 2 . 
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